/* 	ROBUST 2: DEFINE GERMANY WIDE DISPATCH, I.E. NEED ALTERNATIVE FILE FOR LAMBDAS ** 
	* Main adjustments compared to original algorithm: lambda_v1_robust2; ranking
*/   
   

clear all
cd "${PATH}/data"


*******************************
// DEFINE LOCALS // 
local s = $step_size 

* CHOOSE GAMMA // Max share of houses to be covered with solar 
local gamma_min = $gamma_min // 0.05
local gamma_max = $gamma_max // 0.50
local gamma_delta = $gamma_delta // 0.05
local number_gammas = round((`gamma_max' - `gamma_min') / `gamma_delta') + 1

local SCC  $SCC 
*********************************

cap log close 
set logtype text
log using "../results/log_PVoutput_analysis3_10KW_SCC${SCC_n}_robust2.txt", replace


timer clear
timer on 1
set more off

// START DIRECTLY WITH MERGED FILE (see lambdas_v1_robust2)
* Uses common electricity supply accross Germany
use "lambda_15min_robust2.dta", clear



// COMPUTE BASELINE OUTPUT
gen year = yofd(dofc(time_utc))




gen marginal_emissions_EURton = `SCC' * marginal_emissions

// Marginal costs and marginal emissions (mc, me) for each of the steps:
gen mc_step_0 = 0 // Renewables (zero MC)
gen mc_step_1 = 5.1 // Nuclear 
gen mc_step_2 = 6 // Waste and Biomass 
gen mc_step_3 = 11.3  // Lignite
gen mc_step_4 = cost_MWH_coal 
gen mc_step_5 = cost_MWH_ng 
gen mc_step_6 = cost_MWH_crude_oil

gen me_step_0 = 0 // Renewables (zero MC + zero emission)
gen me_step_1 = 0 // Nuclear (no emissions)
gen me_step_2 = 0 // Waste and Biomass (no emissions)
gen me_step_3 = 1.148*`SCC' // Lignite
gen me_step_4 = emission_MWH_coal*`SCC' // hard coal
gen me_step_5 = emission_MWH_ng*`SCC' // gas
gen me_step_6 = emission_MWH_crude_oil*`SCC' // oil

drop cost_* emission_* solar


*Update information on housing stock
gen 	buildings = 3518109 if TSO == "50Hertz"
replace buildings = 6079790 if TSO == "Amprion"
replace buildings = 7034498 if TSO == "TenneT"
replace buildings = 2240361 if TSO == "TransnetBW"

local avg_inst_size = 6.7
gen max_cap = (buildings*`avg_inst_size') /1000 // max cap in MW if all houses had avg. sized solar installation
gen max_cap_gamma = . // max_cap * gamma (in loop)

keep if solar_gen >0 & !missing(solar_gen) // LIMIT sample only to positive solar output hours
****************


** Fraction of small term solar in total production ** 

** 2016 (use max data)
** (FROM solar_capacityTSO.dta, < 10 KW, total cum):
local solarcap_10kw_50Hertz   536.55963
local solarcap_10kw_Amprion   1779.9152
local solarcap_10kw_TenneT    2293.709
local solarcap_10kw_TransnetBW  1117.9264

local fraction_50Hertz   `solarcap_10kw_50Hertz' /   9842.2686
local fraction_Amprion   `solarcap_10kw_Amprion' /   9549.0361
local fraction_TenneT  	 `solarcap_10kw_TenneT'  /   15509.563
local fraction_TransnetBW `solarcap_10kw_TransnetBW' /   5605.6885


// S IS THE SUM OF THE solarcap_10kw_*
local S = 5728 // 5728.1102  // max values of installed capacity (<10KW) in 2016

// Production of residential installations in MWh
rename solar_gen solar_gen_official
gen solar_gen  = .
gen solar_gen_MW = .

levelsof TSO, local(T)
foreach t in `T' {
	replace solar_gen = solar_gen_official * `fraction_`t'' if TSO =="`t'" // production of residential solar systems 
	replace solar_gen_MW = solar_gen_official *`fraction_`t'' /  `solarcap_10kw_`t'' if TSO =="`t'" // production of res. systems [per MW]
}

// For reallocation: solar generation per capacity that is moved around // 
gen solar_gen_s = solar_gen_MW * `s'

*********



// S.D.s FROM LOAD AND SOLAR_SCALED (FROM load_regressions.do)

local factor_sd = 2

local sd_load_50Hertz `factor_sd' * 327.2129
local sd_load_Amprion `factor_sd' * 193.0014
local sd_load_TenneT `factor_sd' * 159.003
local sd_load_TransnetBW `factor_sd' * 93.40471

** Original
local sd_solar_50Hertz `factor_sd' * .0000157
local sd_solar_Amprion `factor_sd' * .0000192
local sd_solar_TenneT `factor_sd' * .000028
local sd_solar_TransnetBW `factor_sd' * .0000219

gen MB_bench  = -dAS_dR + marginal_cost +  marginal_emissions_EURton

*Calculate total value: MB*solar output
gen MB_bench_tot = MB_bench * solar_gen
gen MB_bench_AS = -dAS_dR * solar_gen
gen MB_bench_MC = marginal_cost * solar_gen
gen MB_bench_ME = marginal_emissions_EURton * solar_gen


// AGGREGATE TOTAL BENEFITS 
preserve

	// Collapse by TSO
	collapse (sum) MB_bench_* , by(TSO)
	
	qui: su MB_bench_tot  if TSO == "50Hertz"
		local value_50Hertz_bench = r(mean)
	qui: su  MB_bench_tot   if TSO == "Amprion"	
		local value_Amprion_bench = r(mean)
	qui: su  MB_bench_tot   if TSO == "TenneT"	
		local value_TenneT_bench = r(mean)
	qui: su  MB_bench_tot   if TSO == "TransnetBW"
		local value_TransnetBW_bench = r(mean)
 
	qui: su MB_bench_AS  if TSO == "50Hertz"
		local MB_bench_AS_50Hertz = r(mean)
	qui: su MB_bench_AS  if TSO == "Amprion"
		local MB_bench_AS_Amprion = r(mean)
	qui: su MB_bench_AS  if TSO == "TenneT"
		local MB_bench_AS_TenneT = r(mean)
	qui: su MB_bench_AS  if TSO == "TransnetBW"
		local MB_bench_AS_TransnetBW = r(mean)

	qui: su MB_bench_MC  if TSO == "50Hertz"
		local MB_bench_MC_50Hertz = r(mean)
	qui: su MB_bench_MC  if TSO == "Amprion"
		local MB_bench_MC_Amprion = r(mean)
	qui: su MB_bench_MC  if TSO == "TenneT"
		local MB_bench_MC_TenneT = r(mean)
	qui: su MB_bench_MC  if TSO == "TransnetBW"
		local MB_bench_MC_TransnetBW = r(mean)
		
	qui: su MB_bench_ME  if TSO == "50Hertz"
		local MB_bench_ME_50Hertz = r(mean)
	qui: su MB_bench_ME  if TSO == "Amprion"
		local MB_bench_ME_Amprion = r(mean)
	qui: su MB_bench_ME  if TSO == "TenneT"
		local MB_bench_ME_TenneT = r(mean)
	qui: su MB_bench_ME  if TSO == "TransnetBW"
		local MB_bench_ME_TransnetBW = r(mean)
	
	
	// Collapse for total value
	collapse (sum) MB_bench_* 
	qui: su MB_bench_tot
		local value_bench_total = r(mean)
	qui: su MB_bench_AS
		local value_bench_AS = r(mean)
	qui: su MB_bench_MC
		local value_bench_MC = r(mean)
	qui: su MB_bench_ME
		local value_bench_ME = r(mean)
restore

drop MB_bench MB_bench_* dAS_dR 

*********
// Generate Matrices to store values
*Matrix 1: Allocation, Value. by TSO, total
mat ratios = J(`number_gammas', 20, .)
*Matrix 2: Decomposition of gains: ASC, MC, ME. by TSO, total
mat ratios_decomposition = J(`number_gammas', 32, .)

mat ratios_sd_plus_1_plus_1 = J(`number_gammas', 20, .)
mat ratios_sd_minus_1_minus_1 = J(`number_gammas', 20, .)
mat ratios_sd_minus_1_plus_1 = J(`number_gammas', 20, .)
mat ratios_sd_plus_1_minus_1 = J(`number_gammas', 20, .)



*********
// Reallocation //
cap erase temp_robust2.dta


save temp_robust2.dta

forvalues bands = 0 / 4 {
	
	use temp_robust2, clear
	
	if `bands' == 1 {
		// CASE + +
		local sign_load  1
		local sign_solar 1
	}
	else if `bands' == 2 {
		// CASE - - 
		local sign_load  -1
		local sign_solar -1
	}
	else if `bands' == 3 {
		// CASE + -
		local sign_load  1 
		local sign_solar -1		
	}
	else if `bands' == 4 {
		// CASE - +
		local sign_load  -1
		local sign_solar 1 		
	}

	
	if `bands' > 0 {
		replace load = load + `sign_load' * `sd_load_50Hertz' if TSO == "50Hertz"
		replace load = load + `sign_load' * `sd_load_Amprion' if TSO == "Amprion"
		replace load = load + `sign_load' * `sd_load_TenneT' if TSO == "TenneT"
		replace load = load + `sign_load' * `sd_load_TransnetBW' if TSO == "TransnetBW"
		
		replace solar_gen_MW = cond(solar_gen_MW + `sign_solar' * `sd_solar_50Hertz' > 0, solar_gen_MW + `sign_solar' * `sd_solar_50Hertz', 0) if TSO == "50Hertz"
		replace solar_gen_MW = cond(solar_gen_MW + `sign_solar' * `sd_solar_Amprion' > 0, solar_gen_MW + `sign_solar' * `sd_solar_Amprion', 0) if TSO == "Amprion"
		replace solar_gen_MW = cond(solar_gen_MW + `sign_solar' * `sd_solar_TenneT' > 0, solar_gen_MW + `sign_solar' * `sd_solar_TenneT', 0) if TSO == "TenneT"
		replace solar_gen_MW = cond(solar_gen_MW + `sign_solar' * `sd_solar_TransnetBW' > 0, solar_gen_MW + `sign_solar' * `sd_solar_TransnetBW', 0) if TSO == "TransnetBW"
	
	// NEED TO REDEFINE solar_gen_s
		replace solar_gen_s = solar_gen_MW * `s'
	
	// Check for all load to be positive
		replace load = 0 if load < 0  // in 25 out of 149054 obs. 
	}

	

	local j 1
	* Need to calculate the value for each MW of installed solar and store allocation
	forvalues gamma = `gamma_min'(`gamma_delta')`gamma_max' { 	// LOOP OVER GAMMA VALUES

	
		*Generate variable to keep track of installed in reallocation  
		cap drop solar_cap solar_capl1  cum_solar_gen solar_new value_50Hertz  value_Amprion value_TenneT value_TransnetBW count_excess_load
		cap drop value_AS_* value_MC_* value_ME_*
		gen solar_cap = 0 // solar capacity in loop (i)
		gen solar_capl1 = 0 // solar capacity at iteration (i-1)
		gen cum_solar_gen = 0 // keep track of cumulative solar production -> for correct "step" on supply curve
		gen solar_new = . // define difference of solar_reallocated - solar_original 
		
		gen count_excess_load = 0 // keep track how many times cumulative solar > load
		
		// Store values: Total and decomposition: ASC, MC, ME
		gen value_50Hertz = 0
		gen value_Amprion = 0
		gen value_TenneT = 0
		gen value_TransnetBW = 0
		
		gen value_AS_50Hertz = 0
		gen value_AS_Amprion = 0
		gen value_AS_TenneT = 0
		gen value_AS_TransnetBW = 0
		gen value_MC_50Hertz = 0
		gen value_MC_Amprion = 0
		gen value_MC_TenneT = 0
		gen value_MC_TransnetBW = 0
		gen value_ME_50Hertz = 0
		gen value_ME_Amprion = 0
		gen value_ME_TenneT = 0
		gen value_ME_TransnetBW = 0
		
		// FREQUENCY OF TIMES TO THE RIGHT OR TO THE LEFT OF INITIAL EQM POINT
		local in_flat_MB 0
		local in_step_6 0
		local in_step_5 0
		local in_step_4 0
		local in_step_3 0
		local in_step_2 0
		local in_step_1 0
		local in_step_0 0
		
		
		forvalues i = `s'(`s')`S' {  // loop over all MW of installed capacity to determine optimal allocation
			
			di " " 
			di " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% " 
			di "Reallocating MW: `i' of `S'"
		
			// Set max. capacity per TSO 
			replace max_cap_gamma = max_cap * `gamma'
		
		
			// Cumulative solar generation: include additional MW of solar installed
			replace solar_cap = solar_cap + `s'
			replace cum_solar_gen =  solar_gen_MW * solar_cap
			* solar_gen_s calculated in line 160
		
		
	
*****************************************************
** Consider German-wide bidding zone **
// CALCULATE AS GERMANY WIDE SOLAR PRODUCTION & LOAD (given error bands )
bysort time_utc: egen cum_solar_gen_DE = total(cum_solar_gen)

drop load_DE
bysort time_utc: egen load_DE = total(load)

*****************************************************




// UPDATED MARGINAL EFFECTS ON ANCILLARY SERVICES // 

// Derivatives of above regression model wrt solar
gen dAS_dR = 21.06 +  ( 0.000167) * 2 * cum_solar_gen_DE + ///
	(-4.69e-10) * 3 * cum_solar_gen_DE^2  + (-0.000718) * load_DE + ///
	 ( 5.95e-09) * load_DE^2 + (-2.25e-09) * 2 * cum_solar_gen_DE * load_DE if kload == 1
	

replace dAS_dR = -6.348 + (-0.000231) * 2 * cum_solar_gen_DE + ///
	(-3.21e-09) * 3 * cum_solar_gen_DE^2  + (0.000316) * load_DE + ///
	  (-4.01e-09) * load_DE^2 + (6.37e-09) * 2 * cum_solar_gen_DE * load_DE if kload == 2 		  

replace dAS_dR = 0.273 + 0.0000246  * 2 * cum_solar_gen_DE + ///
	(3.16e-09) * 3 * cum_solar_gen_DE^2  + (-0.0000224) * load_DE + ///
	   3.20e-10  * load_DE^2 + (-1.47e-09) * 2 * cum_solar_gen_DE * load_DE if kload == 3 

				
		
		// Define location on supply curve -> identify MC and marginal Emissions
			replace solar_new =  cum_solar_gen_DE - solar_gen_DE // how much more solar do I generate compared to original allocation? 
			
		
			if (cum_solar_gen_DE <= solar_gen_DE) { 
				*Apply E[MB] as "observed": assumption of "flat MB curves"
				gen MB = -dAS_dR + marginal_cost + marginal_emissions_EURton
				** Track how many times we are in flat part 
				gen MB_tot = cond(cum_solar_gen_DE < load_DE, MB * solar_gen_s, 0) // evaluate MB for block of reallocated solar: 
				gen MB_AS = cond(cum_solar_gen_DE  < load_DE, -dAS_dR * solar_gen_s, 0)
				gen MB_MC = cond(cum_solar_gen_DE  < load_DE, marginal_cost * solar_gen_s, 0)
				gen MB_ME = cond(cum_solar_gen_DE  < load_DE, marginal_emissions_EURton * solar_gen_s , 0)
		
				local in_flat_MB `in_flat_MB' + 1
			}
		
		
			else {
				*We reallocate more solar to TSO as in "benchmark" case -> need to adjust benefits
			
				replace count_excess_load = count_excess_load + 1 // One more interation for which cum solar output > load
			
				if ((solar_new > 0) & (solar_new <= step_6)) {
					* Apply MB from step6
					gen MB = -dAS_dR + mc_step_6 + me_step_6
					
					gen MB_tot = cond(cum_solar_gen_DE  < load_DE, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen_DE  < load_DE, -dAS_dR * solar_gen_s, 0)
					gen MB_MC = cond(cum_solar_gen_DE  < load_DE, mc_step_6 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen_DE  < load_DE, me_step_6 * solar_gen_s, 0)
					local in_step_6 `in_step_6' + 1
		
				}
		
				else if ((solar_new > step_6) & (solar_new <= (step_5 + step_6))) {
					*Apply MB from step 5
					gen MB = -dAS_dR + mc_step_5 + me_step_5
					gen MB_tot = cond(cum_solar_gen_DE  < load_DE, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen_DE  < load_DE, -dAS_dR * solar_gen_s, 0 )
					gen MB_MC = cond(cum_solar_gen_DE  < load_DE, mc_step_5 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen_DE  < load_DE, me_step_5 * solar_gen_s, 0)
					local in_step_5 `in_step_5' + 1
				}
		
				else if  ((solar_new > step_6 + step_5) & (solar_new <= (step_4 + step_5 + step_6))) {
					* Apply MB from step 4
					gen MB = -dAS_dR + mc_step_4 + me_step_4
					gen MB_tot = cond(cum_solar_gen_DE  < load_DE, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen_DE  < load_DE, -dAS_dR * solar_gen_s, 0 )
					gen MB_MC = cond(cum_solar_gen_DE  < load_DE, mc_step_4 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen_DE  < load_DE, me_step_4 * solar_gen_s, 0)
					local in_step_4 `in_step_4' + 1
				}
		
				else if  ((solar_new > step_6 + step_5 + step_4) & (solar_new <= (step_3 + step_4 + step_5 + step_6))) {
					* Apply MB from step 3
					gen MB = -dAS_dR + mc_step_3 + me_step_3
					gen MB_tot = cond(cum_solar_gen_DE  < load_DE, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen_DE  < load_DE, -dAS_dR * solar_gen_s, 0 )
					gen MB_MC = cond(cum_solar_gen_DE  < load_DE, mc_step_3 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen_DE  < load_DE, me_step_3 * solar_gen_s, 0)
					local in_step_3 `in_step_3' + 1
				}
		
				else if  ((solar_new > step_6 + step_5 + step_4 + step_3) & (solar_new <= (step_2 + step_3 + step_4 + step_5 + step_6))) {
					* Apply MB from step 2
					gen MB = -dAS_dR + mc_step_2 + me_step_2
					gen MB_tot = cond(cum_solar_gen_DE  < load_DE, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen_DE  < load_DE, -dAS_dR * solar_gen_s, 0 )
					gen MB_MC = cond(cum_solar_gen_DE  < load_DE, mc_step_2 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen_DE  < load_DE, me_step_2 * solar_gen_s, 0)
					local in_step_2 `in_step_2' + 1
				}
		
				else if  ((solar_new > step_6 + step_5 + step_4 + step_3 + step_2) & (solar_new <= (step_1 + step_2 + step_3 + step_4 + step_5 + step_6))) {
					* Apply MB from step 1
					gen MB = -dAS_dR + mc_step_1 + me_step_1
					gen MB_tot = cond(cum_solar_gen_DE  < load_DE, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen_DE  < load_DE, -dAS_dR * solar_gen_s, 0)
					gen MB_MC = cond(cum_solar_gen_DE  < load_DE, mc_step_1 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen_DE  < load_DE, me_step_1 * solar_gen_s, 0)
					local in_step_1 `in_step_1' + 1
				}
		
				else {
					* Apply MB from step 0
					gen MB = -dAS_dR + mc_step_0 + me_step_0
					gen MB_tot = cond(cum_solar_gen_DE  < load_DE, MB * solar_gen_s, 0) 			
					gen MB_AS = cond(cum_solar_gen_DE  < load_DE, -dAS_dR * solar_gen_s, 0 )
					gen MB_MC = cond(cum_solar_gen_DE  < load_DE, mc_step_0 * solar_gen_s, 0)
					gen MB_ME = cond(cum_solar_gen_DE  < load_DE, me_step_0 * solar_gen_s, 0)
					local in_step_0 `in_step_0' + 1
				}
		
			} //else
		
			
			** OBTAIN RANKING IN AUTOMATED FASHION ** 
			egen meanMB = mean(-MB_tot), by(TSO)
			egen MB_order = axis(meanMB TSO), label(TSO)
			qui: tabstat MB_tot , by(MB_order) s(mean sd N) save
			local first   `r(name1)' 
			local second  `r(name2)' 
			local third   `r(name3)' 
			local fourth  `r(name4)'
			di "RANKING TSO: `first', `second', `third', `fourth'"
	
		
			// Allocation rule: nest options // 
			** Store total solar installed, cumulative value of solar installed
		
			gen alloc_constraint = 1 if (solar_cap > max_cap_gamma) & TSO == "`first'" 
			replace alloc_constraint = 0 if alloc_constraint!=1 // replaces value in all other TSOs
			qui: su alloc_constraint
			local alloc_1 =  `r(mean)'
		
			if `alloc_1' == 0 {
				replace solar_cap = solar_cap - `s' if (TSO == "`second'" | TSO == "`third'" | TSO == "`fourth'")
				* SET back to original value for the TSO's where no allocation - consider constaints!! 
			}
		
		
			else if `alloc_1' != 0 {
		
				drop alloc_constraint
				gen alloc_constraint = 1 if  (solar_cap > max_cap_gamma) & TSO == "`second'" 
				replace alloc_constraint = 0 if alloc_constraint!=1
				qui: su alloc_constraint
				local alloc_2 =  `r(mean)'
		
				if `alloc_2' == 0 {
					replace solar_cap = solar_cap - `s' if (TSO == "`first'" | TSO == "`third'" | TSO == "`fourth'") 
					* SET back to original value for the TSO's where no allocation - consider constaints!! 
				}
		
				else if `alloc_2' != 0 {
		
					drop alloc_constraint
					gen alloc_constraint = 1 if (solar_cap > max_cap_gamma) & TSO == "`third'" 
					replace alloc_constraint = 0 if alloc_constraint!=1
					qui: su alloc_constraint
					local alloc_3 =  `r(mean)'
		
					if `alloc_3' == 0 {
						replace solar_cap = solar_cap - `s' if (TSO == "`first'" | TSO == "`second'" | TSO == "`fourth'") 
						* SET back to original value for the TSO's where no allocation - consider constraints!! 
\					}
		
					else {
						// Put in forth TSO 
						replace solar_cap = solar_cap - `s' if (TSO == "`first'" | TSO == "`second'" | TSO == "`third'") 
					} // TSO 4
				} // TSO 3
			} // TSO 2
		
			// Store total value and in case "decomposition"
			qui: levelsof TSO, local(T)
			foreach t in `T' {
				replace value_`t' = value_`t' + MB_tot if (solar_cap != solar_capl1) & TSO == "`t'"
				replace value_AS_`t' = value_AS_`t' + MB_AS if (solar_cap != solar_capl1) & TSO == "`t'"
				replace value_MC_`t' = value_MC_`t' + MB_MC if (solar_cap != solar_capl1) & TSO == "`t'"
				replace value_ME_`t' = value_ME_`t' + MB_ME if (solar_cap != solar_capl1) & TSO == "`t'"
			}
			// value_TSO: sums total MB for each block of solar that is reallocated
		
			replace solar_capl1 = solar_cap
			drop dAS_dR alloc_constraint MB_order meanMB MB MB_tot MB_AS MB_MC MB_ME
// DROP ALSO CUMULATIVE SOLAR GEN DE // 
			drop cum_solar_gen_DE 

			
		
		} // loop over S (total capacity to reallocate)
		
		
		// AGGREGATE TOTAL BENEFITS 
		qui: su solar_cap if TSO =="50Hertz"
		local solar_cap_50Hertz = r(mean) 
		
		qui: su solar_cap if TSO =="Amprion"
		local solar_cap_Amprion = r(mean)
		
		qui: su solar_cap if TSO =="TenneT"
		local solar_cap_TenneT = r(mean)
		
		qui: su solar_cap if TSO =="TransnetBW"
		local solar_cap_TransnetBW = r(mean)
		
			
		
		// Value by TSO
		
		preserve
		collapse (sum) value* 
		
		qui: su value_50Hertz
			local value_50Hertz = r(mean)
		qui: su value_Amprion
			local value_Amprion = r(mean)
		qui: su value_TenneT
			local value_TenneT = r(mean)
		qui: su  value_TransnetBW
			local value_TransnetBW = r(mean)
			
		// Total of average MB over two year period
		local value_total = `value_50Hertz' + `value_Amprion' + `value_TenneT' + `value_TransnetBW'		
		
		qui: su value_AS_50Hertz  
			local MB_AS_50Hertz = r(mean)
		qui: su value_AS_Amprion  
			local MB_AS_Amprion = r(mean)
		qui: su value_AS_TenneT
			local MB_AS_TenneT = r(mean)
		qui: su value_AS_TransnetBW
			local MB_AS_TransnetBW = r(mean)
		
		qui: su value_MC_50Hertz  
			local MB_MC_50Hertz = r(mean)
		qui: su value_MC_Amprion  
			local MB_MC_Amprion = r(mean)
		qui: su value_MC_TenneT
			local MB_MC_TenneT = r(mean)
		qui: su value_MC_TransnetBW
			local MB_MC_TransnetBW = r(mean)		
		
		qui: su value_ME_50Hertz  
			local MB_ME_50Hertz = r(mean)
		qui: su value_ME_Amprion  
			local MB_ME_Amprion = r(mean)
		qui: su value_ME_TenneT
			local MB_ME_TenneT = r(mean)
		qui: su value_ME_TransnetBW
			local MB_ME_TransnetBW	= r(mean)	
		
		egen value_AS = rowtotal(value_AS_*)
		egen value_MC = rowtotal(value_MC_*)
		egen value_ME = rowtotal(value_ME_*)
		
		qui: su value_AS
			local value_AS_tot = r(mean)
		
		qui: su value_MC
			local value_MC_tot = r(mean)
		
		qui: su value_ME
			local value_ME_tot = r(mean)	
				
				
				
		restore	
			
		di "Bands is `bands'"
		di " "
		
		// Fill in Matrices
	
		if `bands' == 0 {
			* Mat 1: Gains
			mat ratios[`j',1] = `SCC'
			mat ratios[`j',2] = `gamma'
			mat ratios[`j',3] = `solarcap_10kw_50Hertz'
			mat ratios[`j',4] = `solar_cap_50Hertz' 
			mat ratios[`j',5] = `solarcap_10kw_Amprion'
			mat ratios[`j',6] = `solar_cap_Amprion'
			mat ratios[`j',7] = `solarcap_10kw_TenneT'
			mat ratios[`j',8] = `solar_cap_TenneT'
			mat ratios[`j',9] = `solarcap_10kw_TransnetBW'
			mat ratios[`j',10] = `solar_cap_TransnetBW' 
			mat ratios[`j',11] = `value_bench_total' / 1e6
			mat ratios[`j',12] = `value_total' / 1e6
			mat ratios[`j',13] = `value_50Hertz_bench' / 1e6
			mat ratios[`j',14] = `value_50Hertz' / 1e6
			mat ratios[`j',15] = `value_Amprion_bench' / 1e6
			mat ratios[`j',16] = `value_Amprion' / 1e6
			mat ratios[`j',17] = `value_TenneT_bench' / 1e6
			mat ratios[`j',18] = `value_TenneT' / 1e6
			mat ratios[`j',19] = `value_TransnetBW_bench' / 1e6
			mat ratios[`j',20] = `value_TransnetBW' / 1e6
				
			
			*Mat 2: Decomposition
			mat ratios_decomposition[`j',1] = `SCC'
			mat ratios_decomposition[`j',2] = `gamma'
			mat ratios_decomposition[`j',3] = `value_bench_AS' / 1e6
			mat ratios_decomposition[`j',4] = `value_AS_tot' / 1e6
			mat ratios_decomposition[`j',5] = `value_bench_MC' / 1e6
			mat ratios_decomposition[`j',6] =  `value_MC_tot' / 1e6
			mat ratios_decomposition[`j',7] = `value_bench_ME' / 1e6
			mat ratios_decomposition[`j',8] = `value_ME_tot' / 1e6
			mat ratios_decomposition[`j',9] = `MB_bench_AS_50Hertz' / 1e6
			mat ratios_decomposition[`j',10] = `MB_AS_50Hertz' / 1e6
			mat ratios_decomposition[`j',11] = `MB_bench_MC_50Hertz' / 1e6
			mat ratios_decomposition[`j',12] = `MB_MC_50Hertz' / 1e6
			mat ratios_decomposition[`j',13] = `MB_bench_ME_50Hertz' / 1e6
			mat ratios_decomposition[`j',14] = `MB_ME_50Hertz' / 1e6
			mat ratios_decomposition[`j',15] = `MB_bench_AS_Amprion' / 1e6
			mat ratios_decomposition[`j',16] = `MB_AS_Amprion' / 1e6
			mat ratios_decomposition[`j',17] = `MB_bench_MC_Amprion' / 1e6
			mat ratios_decomposition[`j',18] = `MB_MC_Amprion' / 1e6
			mat ratios_decomposition[`j',19] = `MB_bench_ME_Amprion' / 1e6
			mat ratios_decomposition[`j',20] = `MB_ME_Amprion' / 1e6
			mat ratios_decomposition[`j',21] = `MB_bench_AS_TenneT' / 1e6
			mat ratios_decomposition[`j',22] = `MB_AS_TenneT' / 1e6
			mat ratios_decomposition[`j',23] = `MB_bench_MC_TenneT' / 1e6
			mat ratios_decomposition[`j',24] = `MB_MC_TenneT' / 1e6
			mat ratios_decomposition[`j',25] = `MB_bench_ME_TenneT' / 1e6
			mat ratios_decomposition[`j',26] = `MB_ME_TenneT' / 1e6
			mat ratios_decomposition[`j',27] = `MB_bench_AS_TransnetBW' / 1e6
			mat ratios_decomposition[`j',28] = `MB_AS_TransnetBW' / 1e6
			mat ratios_decomposition[`j',29] = `MB_bench_MC_TransnetBW' / 1e6
			mat ratios_decomposition[`j',30] = `MB_MC_TransnetBW' / 1e6
			mat ratios_decomposition[`j',31] = `MB_bench_ME_TransnetBW' / 1e6
			mat ratios_decomposition[`j',32] = `MB_ME_TransnetBW' / 1e6
			** ALL VALUES (value, MB, expressed in Million EUR)
		}	
		else if `bands' == 1 {
			* Mat 3: Gains 
			mat ratios_sd_plus_1_plus_1[`j',1] = `SCC'
			mat ratios_sd_plus_1_plus_1[`j',2] = `gamma'
			mat ratios_sd_plus_1_plus_1[`j',3] = `solarcap_10kw_50Hertz'
			mat ratios_sd_plus_1_plus_1[`j',4] = `solar_cap_50Hertz' 
			mat ratios_sd_plus_1_plus_1[`j',5] = `solarcap_10kw_Amprion'
			mat ratios_sd_plus_1_plus_1[`j',6] = `solar_cap_Amprion'
			mat ratios_sd_plus_1_plus_1[`j',7] = `solarcap_10kw_TenneT'
			mat ratios_sd_plus_1_plus_1[`j',8] = `solar_cap_TenneT'
			mat ratios_sd_plus_1_plus_1[`j',9] = `solarcap_10kw_TransnetBW'
			mat ratios_sd_plus_1_plus_1[`j',10] = `solar_cap_TransnetBW' 
			mat ratios_sd_plus_1_plus_1[`j',11] = `value_bench_total' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',12] = `value_total' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',13] = `value_50Hertz_bench' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',14] = `value_50Hertz' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',15] = `value_Amprion_bench' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',16] = `value_Amprion' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',17] = `value_TenneT_bench' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',18] = `value_TenneT' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',19] = `value_TransnetBW_bench' / 1e6
			mat ratios_sd_plus_1_plus_1[`j',20] = `value_TransnetBW' / 1e6
		}
		else if `bands' == 2 {
			* Mat 4: Gains 
			mat ratios_sd_minus_1_minus_1[`j',1] = `SCC'
			mat ratios_sd_minus_1_minus_1[`j',2] = `gamma'
			mat ratios_sd_minus_1_minus_1[`j',3] = `solarcap_10kw_50Hertz'
			mat ratios_sd_minus_1_minus_1[`j',4] = `solar_cap_50Hertz' 
			mat ratios_sd_minus_1_minus_1[`j',5] = `solarcap_10kw_Amprion'
			mat ratios_sd_minus_1_minus_1[`j',6] = `solar_cap_Amprion'
			mat ratios_sd_minus_1_minus_1[`j',7] = `solarcap_10kw_TenneT'
			mat ratios_sd_minus_1_minus_1[`j',8] = `solar_cap_TenneT'
			mat ratios_sd_minus_1_minus_1[`j',9] = `solarcap_10kw_TransnetBW'
			mat ratios_sd_minus_1_minus_1[`j',10] = `solar_cap_TransnetBW' 
			mat ratios_sd_minus_1_minus_1[`j',11] = `value_bench_total' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',12] = `value_total' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',13] = `value_50Hertz_bench' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',14] = `value_50Hertz' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',15] = `value_Amprion_bench' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',16] = `value_Amprion' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',17] = `value_TenneT_bench' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',18] = `value_TenneT' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',19] = `value_TransnetBW_bench' / 1e6
			mat ratios_sd_minus_1_minus_1[`j',20] = `value_TransnetBW' / 1e6
		}
		else if `bands' == 3 {
			* Mat 4: Gains 
			mat ratios_sd_plus_1_minus_1[`j',1] = `SCC'
			mat ratios_sd_plus_1_minus_1[`j',2] = `gamma'
			mat ratios_sd_plus_1_minus_1[`j',3] = `solarcap_10kw_50Hertz'
			mat ratios_sd_plus_1_minus_1[`j',4] = `solar_cap_50Hertz' 
			mat ratios_sd_plus_1_minus_1[`j',5] = `solarcap_10kw_Amprion'
			mat ratios_sd_plus_1_minus_1[`j',6] = `solar_cap_Amprion'
			mat ratios_sd_plus_1_minus_1[`j',7] = `solarcap_10kw_TenneT'
			mat ratios_sd_plus_1_minus_1[`j',8] = `solar_cap_TenneT'
			mat ratios_sd_plus_1_minus_1[`j',9] = `solarcap_10kw_TransnetBW'
			mat ratios_sd_plus_1_minus_1[`j',10] = `solar_cap_TransnetBW' 
			mat ratios_sd_plus_1_minus_1[`j',11] = `value_bench_total' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',12] = `value_total' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',13] = `value_50Hertz_bench' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',14] = `value_50Hertz' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',15] = `value_Amprion_bench' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',16] = `value_Amprion' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',17] = `value_TenneT_bench' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',18] = `value_TenneT' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',19] = `value_TransnetBW_bench' / 1e6
			mat ratios_sd_plus_1_minus_1[`j',20] = `value_TransnetBW' / 1e6
		}
		else if `bands' == 4 {
			* Mat 4: Gains 
			mat ratios_sd_minus_1_plus_1[`j',1] = `SCC'
			mat ratios_sd_minus_1_plus_1[`j',2] = `gamma'
			mat ratios_sd_minus_1_plus_1[`j',3] = `solarcap_10kw_50Hertz'
			mat ratios_sd_minus_1_plus_1[`j',4] = `solar_cap_50Hertz' 
			mat ratios_sd_minus_1_plus_1[`j',5] = `solarcap_10kw_Amprion'
			mat ratios_sd_minus_1_plus_1[`j',6] = `solar_cap_Amprion'
			mat ratios_sd_minus_1_plus_1[`j',7] = `solarcap_10kw_TenneT'
			mat ratios_sd_minus_1_plus_1[`j',8] = `solar_cap_TenneT'
			mat ratios_sd_minus_1_plus_1[`j',9] = `solarcap_10kw_TransnetBW'
			mat ratios_sd_minus_1_plus_1[`j',10] = `solar_cap_TransnetBW' 
			mat ratios_sd_minus_1_plus_1[`j',11] = `value_bench_total' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',12] = `value_total' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',13] = `value_50Hertz_bench' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',14] = `value_50Hertz' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',15] = `value_Amprion_bench' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',16] = `value_Amprion' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',17] = `value_TenneT_bench' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',18] = `value_TenneT' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',19] = `value_TransnetBW_bench' / 1e6
			mat ratios_sd_minus_1_plus_1[`j',20] = `value_TransnetBW' / 1e6
		}
				
		// CHECK LOCATION ON SUPPLY CURVE
		di " "
		di " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "
		di " "
		di "FREQUENCY OF TIMES TO THE LEFT AND TO THE RIGHT OF INITIAL EQM POINT"
		di `in_flat_MB'
		di `in_step_6'
		di `in_step_5'
		di `in_step_4'
		di `in_step_3'
		di `in_step_2'
		di `in_step_1'
		di `in_step_0'
			
		di " "
		di " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "
		di "Number of times solar production is larger than load "
		qui: su count_excess_load if TSO == "50Hertz"
		di "50 Hertz: " `r(mean)'
		qui: su count_excess_load if TSO == "TenneT"
		di "TenneT: " `r(mean)'
		qui: su count_excess_load if TSO == "Amprion"
		di "Amprion: " `r(mean)'
		qui: su count_excess_load if TSO == "TransnetBW"
		di "TransnetBW: " `r(mean)'
		
		
		local ++j
		di `j'

	} // LOOP OVER GAMMA

	
} // LOOP OVER BANDS	




// Extract values + plot
mat colnames ratios = SCC gamma cap_50HZ cap_50HZ_alloc cap_Amprion cap_Amprion_alloc cap_TenneT cap_TenneT_alloc cap_Transnet cap_Transnet_alloc v_bench v_alloc /// 
v_50HZ v_50HZ_alloc v_Amprion v_Amprion_alloc v_TenneT v_TenneT_alloc v_Transnet v_Transnet_alloc
svmat ratios, names(col)
	
mat colnames ratios_decomposition = SCC_ gamma_ v_AS v_AS_alloc v_MC v_MC_alloc v_ME v_ME_alloc v_AS_50HZ v_AS_50HZ_alloc v_MC_50HZ v_MC_50HZ_alloc ///
v_ME_50HZ v_ME_50HZ_alloc v_AS_Amp v_AS_50_HZ_alloc v_MC_Amp v_MC_Amp_alloc v_ME_Amp v_ME_Amp_alloc ///
v_AS_TenneT v_AS_TenneT_alloc v_MC_TenneT v_MC_TenneT_alloc v_ME_TenneT v_ME_TenneT_alloc v_AS_Trans v_AS_Trans_alloc v_MC_Trans v_MC_Trans_alloc v_ME_Trans v_ME_Trans_alloc
svmat ratios_decomposition, names(col)
	
*
mat colnames ratios_sd_plus_1_plus_1 = SCC_sd_p_1_p_1 gamma_sd_p_1_p_1 cap_50HZ_sd_p_1_p_1 cap_50HZ_alloc_sd_p_1_p_1 cap_Amprion_sd_p_1_p_1 cap_Amprion_alloc_sd_p_1_p_1 cap_TenneT_sd_p_1_p_1 cap_TenneT_alloc_sd_p_1_p_1 cap_Transnet_sd_p_1_p_1 cap_Transnet_alloc_sd_p_1_p_1 v_bench_sd_p_1_p_1 v_alloc_sd_p_1_p_1 /// 
v_50HZ_sd_p_1_p_1 v_50HZ_alloc_sd_p_1_p_1 v_Amprion_sd_p_1_p_1 v_Amprion_alloc_sd_p_1_p_1 v_TenneT_sd_p_1_p_1 v_TenneT_alloc_sd_p_1_p_1 v_Transnet_sd_p_1_p_1 v_Transnet_alloc_sd_p_1_p_1
svmat ratios_sd_plus_1_plus_1, names(col)		

mat colnames ratios_sd_minus_1_minus_1 = SCC_sd_m_1_m_1 gamma_sd_m_1_m_1 cap_50HZ_sd_m_1_m_1 cap_50HZ_alloc_sd_m_1_m_1 cap_Amprion_sd_m_1_m_1 cap_Amprion_alloc_sd_m_1_m_1 cap_TenneT_sd_m_1_m_1 cap_TenneT_alloc_sd_m_1_m_1 cap_Transnet_sd_m_1_m_1 cap_Transnet_alloc_sd_m_1_m_1 v_bench_sd_m_1_m_1 v_alloc_sd_m_1_m_1 /// 
v_50HZ_sd_m_1_m_1 v_50HZ_alloc_sd_m_1_m_1 v_Amprion_sd_m_1_m_1 v_Amprion_alloc_sd_m_1_m_1 v_TenneT_sd_m_1_m_1 v_TenneT_alloc_sd_m_1_m_1 v_Transnet_sd_m_1_m_1 v_Transnet_alloc_sd_m_1_m_1
svmat ratios_sd_minus_1_minus_1, names(col)		

mat colnames ratios_sd_plus_1_minus_1 = SCC_sd_p_1_m_1 gamma_sd_p_1_m_1 cap_50HZ_sd_p_1_m_1 cap_50HZ_alloc_sd_p_1_m_1 cap_Amprion_sd_p_1_m_1 cap_Amprion_alloc_sd_p_1_m_1 cap_TenneT_sd_p_1_m_1 cap_TenneT_alloc_sd_p_1_m_1 cap_Transnet_sd_p_1_m_1 cap_Transnet_alloc_sd_p_1_m_1 v_bench_sd_p_1_m_1 v_alloc_sd_p_1_m_1 /// 
v_50HZ_sd_p_1_m_1 v_50HZ_alloc_sd_p_1_m_1 v_Amprion_sd_p_1_m_1 v_Amprion_alloc_sd_p_1_m_1 v_TenneT_sd_p_1_m_1 v_TenneT_alloc_sd_p_1_m_1 v_Transnet_sd_p_1_m_1 v_Transnet_alloc_sd_p_1_m_1
svmat ratios_sd_plus_1_minus_1, names(col)		

mat colnames ratios_sd_minus_1_plus_1 = SCC_sd_m_1_p_1 gamma_sd_m_1_p_1 cap_50HZ_sd_m_1_p_1 cap_50HZ_alloc_sd_m_1_p_1 cap_Amprion_sd_m_1_p_1 cap_Amprion_alloc_sd_m_1_p_1 cap_TenneT_sd_m_1_p_1 cap_TenneT_alloc_sd_m_1_p_1 cap_Transnet_sd_m_1_p_1 cap_Transnet_alloc_sd_m_1_p_1 v_bench_sd_m_1_p_1 v_alloc_sd_m_1_p_1 /// 
v_50HZ_sd_m_1_p_1 v_50HZ_alloc_sd_m_1_p_1 v_Amprion_sd_m_1_p_1 v_Amprion_alloc_sd_m_1_p_1 v_TenneT_sd_m_1_p_1 v_TenneT_alloc_sd_m_1_p_1 v_Transnet_sd_m_1_p_1 v_Transnet_alloc_sd_m_1_p_1
svmat ratios_sd_minus_1_plus_1, names(col)		
*

log close 

// Erase temp file
cap erase temp_robust2.dta

// STORE MAIN VALUES
preserve
keep SCC- v_Transnet_alloc_sd_m_1_p_1 
drop if missing(SCC)
save "PVoutput_analysis3_10KW_v3_SCC${SCC_n}_robust2.dta", replace
restore 
